suppressWarnings(library(mvtnorm))
suppressWarnings(library(R2HTML))
suppressWarnings(library(npmc))
# retrieve args
Args <- commandArgs(TRUE)

#Read in data
statdata <- read.csv(Args[3], header=TRUE, sep=",")

#Copy Args
response <- Args[4]
treatment <- Args[5]
sig <- 1 - as.numeric(Args[6])
statstest <- Args[7]
contlevel <- Args[8]

#Create variable to use as a formula
respVTreat <- paste(response, "~", treatment, sep="")

#Setup the html file and associated css file
htmlFile <- sub(".csv", ".html", Args[3]); #determine the file name of the html file
.HTML.file = htmlFile
cssFile <- "r2html.css"
cssFile <- paste("'",cssFile,"'", sep="") #need to enclose in quotes when path has spaces in it
HTMLCSS(CSSfile = cssFile)

#Output HTML header
HTML.title("<bf>SilveR Non-Parametric Analysis", HR=1, align="left")

#Removing illegal charaters

YAxisTitle<-response
XAxisTitle<-treatment

for (i in 1:10)
{

# Additional characters included Aug 2010 (STB)
YAxisTitle<-sub("ivs_tilde_ivs"	,"~", YAxisTitle) 
YAxisTitle<-sub("ivs_star_ivs"	,"*", YAxisTitle) 
YAxisTitle<-sub("ivs_plus_ivs"	,"+", YAxisTitle) 

YAxisTitle<-sub("ivs_sp_ivs"	," ", YAxisTitle) 
YAxisTitle<-sub("ivs_ob_ivs"	,"(", YAxisTitle) 
YAxisTitle<-sub("ivs_cb_ivs"	,")", YAxisTitle) 
YAxisTitle<-sub("ivs_div_ivs"	,"/", YAxisTitle) 
YAxisTitle<-sub("ivs_pc_ivs"	,"%", YAxisTitle) 
YAxisTitle<-sub("ivs_hash_ivs"	,"#", YAxisTitle) 
YAxisTitle<-sub("ivs_pt_ivs"	,".", YAxisTitle) 
YAxisTitle<-sub("ivs_hyphen_ivs","-", YAxisTitle) 
YAxisTitle<-sub("ivs_at_ivs"	,"@", YAxisTitle) 
YAxisTitle<-sub("ivs_colon_ivs"	,":", YAxisTitle) 
YAxisTitle<-sub("ivs_exclam_ivs","!", YAxisTitle) 
YAxisTitle<-sub("ivs_quote_ivs"	,"`", YAxisTitle) 
YAxisTitle<-sub("ivs_pound_ivs"	,"£", YAxisTitle) 
YAxisTitle<-sub("ivs_dollar_ivs","$", YAxisTitle) 
YAxisTitle<-sub("ivs_hat_ivs"	,"^", YAxisTitle) 
YAxisTitle<-sub("ivs_amper_ivs"	,"&", YAxisTitle) 
YAxisTitle<-sub("ivs_obrace_ivs","{", YAxisTitle) 
YAxisTitle<-sub("ivs_cbrace_ivs","}", YAxisTitle) 
YAxisTitle<-sub("ivs_semi_ivs"	,";", YAxisTitle) 
YAxisTitle<-sub("ivs_pipe_ivs"	,"|", YAxisTitle) 
YAxisTitle<-sub("ivs_slash_ivs"	,"\\", YAxisTitle) 
YAxisTitle<-sub("ivs_osb_ivs"	,"[", YAxisTitle) 
YAxisTitle<-sub("ivs_csb_ivs"	,"]", YAxisTitle) 
YAxisTitle<-sub("ivs_eq_ivs"	,"=", YAxisTitle) 
YAxisTitle<-sub("ivs_lt_ivs"	,"<", YAxisTitle) 
YAxisTitle<-sub("ivs_gt_ivs"	,">", YAxisTitle) 
YAxisTitle<-sub("ivs_dblquote_ivs"	,"\"", YAxisTitle) 


# Additional characters included Aug 2010 (STB)
XAxisTitle<-sub("ivs_tilde_ivs"	,"~", XAxisTitle) 
XAxisTitle<-sub("ivs_star_ivs"	,"*", XAxisTitle) 
XAxisTitle<-sub("ivs_plus_ivs"	,"+", XAxisTitle) 

XAxisTitle<-sub("ivs_sp_ivs"	," ", XAxisTitle) 
XAxisTitle<-sub("ivs_ob_ivs"	,"(", XAxisTitle) 
XAxisTitle<-sub("ivs_cb_ivs"	,")", XAxisTitle) 
XAxisTitle<-sub("ivs_div_ivs"	,"/", XAxisTitle) 
XAxisTitle<-sub("ivs_pc_ivs"	,"%", XAxisTitle) 
XAxisTitle<-sub("ivs_hash_ivs"	,"#", XAxisTitle) 
XAxisTitle<-sub("ivs_pt_ivs"	,".", XAxisTitle) 
XAxisTitle<-sub("ivs_hyphen_ivs","-", XAxisTitle) 
XAxisTitle<-sub("ivs_at_ivs"	,"@", XAxisTitle) 
XAxisTitle<-sub("ivs_colon_ivs"	,":", XAxisTitle) 
XAxisTitle<-sub("ivs_exclam_ivs","!", XAxisTitle) 
XAxisTitle<-sub("ivs_quote_ivs"	,"`", XAxisTitle) 
XAxisTitle<-sub("ivs_pound_ivs"	,"£", XAxisTitle) 
XAxisTitle<-sub("ivs_dollar_ivs","$", XAxisTitle) 
XAxisTitle<-sub("ivs_hat_ivs"	,"^", XAxisTitle) 
XAxisTitle<-sub("ivs_amper_ivs"	,"&", XAxisTitle) 
XAxisTitle<-sub("ivs_obrace_ivs","{", XAxisTitle) 
XAxisTitle<-sub("ivs_cbrace_ivs","}", XAxisTitle) 
XAxisTitle<-sub("ivs_semi_ivs"	,";", XAxisTitle) 
XAxisTitle<-sub("ivs_pipe_ivs"	,"|", XAxisTitle) 
XAxisTitle<-sub("ivs_slash_ivs"	,"\\", XAxisTitle) 
XAxisTitle<-sub("ivs_osb_ivs"	,"[", XAxisTitle) 
XAxisTitle<-sub("ivs_csb_ivs"	,"]", XAxisTitle) 
XAxisTitle<-sub("ivs_eq_ivs"	,"=", XAxisTitle) 
XAxisTitle<-sub("ivs_lt_ivs"	,"<", XAxisTitle) 
XAxisTitle<-sub("ivs_gt_ivs"	,">", XAxisTitle) 
XAxisTitle<-sub("ivs_dblquote_ivs"	,"\"", XAxisTitle) 
}

#Response

title<-c("Response")
HTML.title(title, HR=2, align="left")
add<-paste(c("The  "), response, sep="")
add<-paste(add, " response is currently being analysed by the Non-Parametric Analysis module", sep="")

if(treatment != "NULL")
{
	add<-paste(add, c(", with  "), sep="")
	add<-paste(add, treatment, sep="")
	add<-paste(add, " fitted as the treatment factor.", sep="")
} else {
	add<-paste(add, ".", sep="")
}

HTML.title("</bf> ", HR=2, align="left")
HTML.title(add, HR=0, align="left")






#set-up blank table entry
blank<-c(" ")
index<-1

#Pull out the columns we are interested in, get rid of missing data, and convert the treatment to character level (if numeric)
temp <- list(factor(as.character(eval(parse(text = paste("statdata$", treatment))))), eval(parse(text = paste("statdata$", response))))
statdata <- data.frame(temp)

cols<-c(treatment, response)
colnames(statdata)<-cols
statdata<-na.omit(statdata)
leng<-length(levels(as.factor(eval(parse(text = paste("statdata$", treatment))))))

#Basic stats by group
Minimum<-c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
Maximum<-c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
Q1<-c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
Median<-c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
Q3<-c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
 
length(Minimum)=length(unique(eval(parse(text = paste("statdata$", treatment)))))
length(Maximum)=length(unique(eval(parse(text = paste("statdata$", treatment)))))
length(Q1)=length(unique(eval(parse(text = paste("statdata$", treatment)))))
length(Median)=length(unique(eval(parse(text = paste("statdata$", treatment)))))
length(Q3)=length(unique(eval(parse(text = paste("statdata$", treatment)))))
	
for(i in 1:length(unique(eval(parse(text = paste("statdata$", treatment)))))){Minimum[i]<-min(split(eval(parse(text = paste("statdata$", response))),eval(parse(text = paste("statdata$", treatment))))[[i]])}
for(i in 1:length(unique(eval(parse(text = paste("statdata$", treatment)))))){Q1[i]<-boxplot.stats(split(eval(parse(text = paste("statdata$", response))),eval(parse(text = paste("statdata$", treatment))))[[i]])$stats[2]}
for(i in 1:length(unique(eval(parse(text = paste("statdata$", treatment)))))){Median[i]<-boxplot.stats(split(eval(parse(text = paste("statdata$", response))),eval(parse(text = paste("statdata$", treatment))))[[i]])$stats[3]}
for(i in 1:length(unique(eval(parse(text = paste("statdata$", treatment)))))){Q3[i]<-boxplot.stats(split(eval(parse(text = paste("statdata$", response))),eval(parse(text = paste("statdata$", treatment))))[[i]])$stats[4]}
for(i in 1:length(unique(eval(parse(text = paste("statdata$", treatment)))))){Maximum[i]<-max(split(eval(parse(text = paste("statdata$", response))),eval(parse(text = paste("statdata$", treatment))))[[i]])}

resultz<-cbind(Minimum,Q1, Median, Q3,Maximum)
results<-format(round(resultz,2),nsmall=2)

header<-c(" ", " "," " , " ", " ")
tables<-rbind(header, results)

good1<-c("Group")

fort<-levels(eval(parse(text = paste("statdata$", treatment))))
for(i in 1:length(unique(eval(parse(text = paste("statdata$", treatment)))))){good1[i+1]<-fort[i]}

row.names(tables)<-good1

table<-data.frame(tables)
HTMLbr()
HTML.title("<bf>Summary data", HR=2, align="left")
HTML(table, classfirstline="second", align="left")

#Boxplot of non-para data
title<-paste("Box-plot of ", response, sep="")
title<-paste(title, ", categorised by ", sep="")
title<-paste(title, treatment, sep="")

plotFile <- sub(".html", "boxplot.jpg", htmlFile)
jpeg(plotFile)


#boxplot(as.formula(eval(parse(text = paste("statdata$", response)))~factor(eval(parse(text = paste("statdata$", treatment))))),  main=title)
bw<-	boxplot(as.formula(eval(parse(text = paste("statdata$", response)))~factor(eval(parse(text = paste("statdata$", treatment))))),  main=title, col="black", range = 1.5, outline=TRUE)
par(las=1)
bxp(bw,boxfill= "grey95", whisklty=1,medpch=16, medlty="blank",xlab=XAxisTitle, 	ylab =YAxisTitle)

HTMLbr()
HTML.title("<bf>Box-plot", HR=2, align="left")
void<-HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", plotFile), Align="left")

#Main Analysis
if (statstest=="MannWhitney" || (statstest=="AllComparisons" && leng == 2 ) || (statstest=="CompareToControl" && leng == 2 ))
{
	#Kruskal Wallis and Mann Whitney
	dim <- length(unique(eval(parse(text = paste("statdata$", treatment)))))
	
	if (dim > 2)
	{
		HTMLbr()
		HTML.title("<bf>Kruskal-Wallis test", HR=2, align="left")
		options<-(scipen=20)
		kruskalOut <- kruskal.test(as.formula(respVTreat), data=statdata)

		pvalue<-format(round(kruskalOut$p.value, 4), nsmall=4, scientific=FALSE)

		for (i in 1:(length(pvalue))) 
		{
			if (pvalue[i]<0.001) 
			{
				pvalue[i]<-0.001
				pvalue[i]<- paste("<",pvalue[i])
			}
		}

		Wstat<- format(round(kruskalOut$statistic,2),nsmall=2)
		df<- format(round(kruskalOut$parameter,0),nsmall=0)
		colname<-c("Test statistic", "DF", "p-value")
		rowname<-c("Result")
		temptab<-cbind(Wstat, df, pvalue)
		temptable<-data.frame(temptab)
		colnames(temptable)<-colname
		row.names(temptable)<-rowname
		HTML(temptable, classfirstline="second", align="left")				
	} else if (dim==2)
	{
		HTMLbr()
		HTML.title("<bf> Mann-Whitney test", HR=2, align="left")		
		mess<-c("You have selected a factor with only two levels, hence a Mann-Whitney test, also know as a Wilcoxon rank sum test, has been used to analyse the data rather than a Kruskal Wallis test.")
		HTML.title("</bf> ", HR=2, align="left")
		HTML.title(mess, HR=0, align="left")

	
		wilcoxOut <- wilcox.test(as.formula(respVTreat), data=statdata, exact=TRUE, conf.int = TRUE, conf.level = sig)
		pvalue<-format(round(wilcoxOut$p.value, 4), nsmall=4, scientific=FALSE)

		for (i in 1:(length(pvalue))) 
		{
			if (pvalue[i]<0.001) 
			{
				pvalue[i]<-0.001
				pvalue[i]<- paste("<",pvalue[i])
			}
		}

		Wstat<- format(round(wilcoxOut$statistic,2),nsmall=2)
		temptab<-cbind(Wstat, pvalue)
		temptable<-data.frame(temptab)
		colname<-c("Test statistic", "p-value")
		colnames(temptable)<-colname
		rowname<-c("Test result")
		row.names(temptable)<-rowname
		HTML(temptable, classfirstline="second", align="left")

		if(length(unique(eval(parse(text = paste("statdata$", response)))))==length((eval(parse(text = paste("statdata$", treatment))))))
		{
			HTML.title("</bf> ", HR=2, align="left")
			HTML.title("<bf>As there are no ties in the responses, the exact test result has been calculated.", HR=0, align="left")
		}
		if(length(unique(eval(parse(text = paste("statdata$", response))))) < length((eval(parse(text = paste("statdata$", treatment))))))
		{
			HTML.title("</bf> ", HR=2, align="left")
			HTML.title("<bf>As there are ties in the responses, the asymptotic test result has been calculated.", HR=0, align="left")
		}
	}
} else if (statstest=="AllComparisons" && leng >= 3)
{
	#All comparisons test
	npmcallcomp<-function(response, treatment){
	dataset<-data.frame(var=response, class=treatment)
	summary(npmc(dataset, df=1, alpha=sig),type="BF")}

	allpair<-npmcallcomp(eval(parse(text = paste("statdata$", response))), eval(parse(text = paste("statdata$", treatment))))
    	Treatment<-allpair$'Data-structure'[2]
	Nobs<-allpair$'Data-structure'[3]

	temp3<-cbind(Treatment, Nobs)
	header<-c(NA, " ")
	allpairtable<-rbind(header, temp3)
	allpairtable1<-data.frame(allpairtable)

	temp4<-c("Treatment", "Nobs")
	colnames(allpairtable1)<-temp4

	good2<-c("Treatment ID")
	for(i in 1:length(unique(eval(parse(text = paste("statdata$", treatment)))))){good2[i+1]<-i}
	row.names(allpairtable1)<-good2

	# Behrens Fisher results
	HTMLbr()
	HTML.title("<bf>All pairwise comparisons - Behrens Fisher tests", HR=2, align="left")

	compno<-c(0,1,3,6,10,15,21,28,36,45,55,66,78,91,105,120,136,153,171,190,210,231,253,276,300,325,351,378,406,435)
	int<-compno[length(unique(eval(parse(text = paste("statdata$", treatment)))))]+1

	Comparison<-allpair$'Results of the multiple Behrens-Fisher-Test'[1]
	Pvalue<-format(round(allpair$'Results of the multiple Behrens-Fisher-Test'[6],4),nsmall=4, scientific=FALSE)

	y<-1
	allgroup1<-unique(eval(parse(text = paste("statdata$", treatment))))
	for(i in 1:(length(unique(eval(parse(text = paste("statdata$", treatment)))))-1)) 
	{
		for(j in 1:(length(unique(eval(parse(text = paste("statdata$", treatment)))))-i)) 
		{
			allgroup1[y]<-unique(eval(parse(text = paste("statdata$", treatment))))[i]
			y<-y+1
		}
	}

	allgroup2<-unique(eval(parse(text = paste("statdata$", treatment))))
	z<-1
	for(i in 1:(length(unique(eval(parse(text = paste("statdata$", treatment)))))-1)) 
	{
		for(j in 1:(length(unique(eval(parse(text = paste("statdata$", treatment)))))-i)) 
		{
			allgroup2[z]<-unique(eval(parse(text = paste("statdata$", treatment))))[j+i]
			z<-z+1
		}
	}
	allvs<-c("r")	
	for(i in 1:(int-1)){allvs[i]<-"vs."}

	temp5<-cbind(allgroup1, allvs, allgroup2, Pvalue)

	for (i in 1:(dim(temp5)[1])) 
	{
		if (temp5[i,4]<0.001) 
		{
			temp5[i,4]<-0.001
			temp5[i,4]<- paste("<",temp5[i,4])
		}
	}

	header<-c(NA, NA, NA, NA)

	allpairtab<-rbind(header, temp5)
	allpairtable<-data.frame(allpairtab)

	temp6<-c("Gp 1", "vs.", "Gp 2", "p-value")
	colnames(allpairtable)<-temp6

	good3<-c("Comparison")
	for(i in 2:int){good3[i]<-i-1}
	row.names(allpairtable)<-good3
	
	HTML(allpairtable, classfirstline="second", align="left")

	compno<-c(0,1,3,6,10,15,21,28,36,45,55,66,78,91,105,120,136,153,171,190,210,231,253,276,300,325,351,378,406,435)
	int<-compno[length(unique(eval(parse(text = paste("statdata$", treatment)))))]+1
	allvs<-c("r")	
	for(i in 1:(int-1)){allvs[i]<-"vs."}

	y<-1
	allgroup1<-unique(eval(parse(text = paste("statdata$", treatment))))
	for(i in 1:(length(unique(eval(parse(text = paste("statdata$", treatment)))))-1)) 
	{
		for(j in 1:(length(unique(eval(parse(text = paste("statdata$", treatment)))))-i)) 
		{
			allgroup1[y]<-unique(eval(parse(text = paste("statdata$", treatment))))[i]
			y<-y+1
		}
	}

	allgroup2<-unique(eval(parse(text = paste("statdata$", treatment))))
	z<-1
	for(i in 1:(length(unique(eval(parse(text = paste("statdata$", treatment)))))-1)) 
	{
		for(j in 1:(length(unique(eval(parse(text = paste("statdata$", treatment)))))-i)) 
		{
			allgroup2[z]<-unique(eval(parse(text = paste("statdata$", treatment))))[j+i]
			z<-z+1
		}
	}
	
	tabletemp<-data.frame(allgroup1, allvs, allgroup2)
	blankz<-c(NA,NA,NA)
	tabletemp<-rbind(blankz,tabletemp)	

	for(s in 1:(length(unique(eval(parse(text = paste("statdata$", treatment)))))-1)) 
	{
		for(t in s+1:(length(unique(eval(parse(text = paste("statdata$", treatment)))))-s)) 
		{
			# seperate out groups
			data1<-c(1)
			data1<-split(eval(parse(text = paste("statdata$",response))),eval(parse(text = paste("statdata$", treatment))))[[s]]
			data2<-c(1)
			data2<-split(eval(parse(text = paste("statdata$",response))),eval(parse(text = paste("statdata$", treatment))))[[t]]
			resp<-c(data1,data2)
	
			#Set up the factor names column
			length1<-length(data1)
			length2<-length(data2)
			fact1<-c("a")
			fact2<-c("b")
			for(m in 1:length1) {fact1[m]<- unique(eval(parse(text = paste("statdata$", treatment))))[[s]]}
			for(n in 1:length2) {fact2[n]<- unique(eval(parse(text = paste("statdata$", treatment))))[[t]]}
			grou<-c(fact1,fact2)
	
			# combine the datasets
			comb1<-cbind(grou,resp)
			finaldata<-data.frame(comb1)

			#Wilcoxon test
			wilcox<-wilcox.test(resp~grou, exact=FALSE)
			pv<-format(round(wilcox$p.value, 4), nsmall=4, scientific=FALSE)

			for (i in 1:(length(pv))) 
			{
				if (pv[i]<0.001) 
				{
					pv[i]<-0.001
					pv[i]<- paste("<",pv[i])
				}
			}

			#merge results
			lines<-c(pv)
			blank<-rbind(blank, lines)
	
			index<-index+1
		}
	}
		
	blank2<-blank
	blank<-cbind(tabletemp,blank)
		
	#create final table
	temp16<-c("Gp 1", "vs.", "Gp 2", "p-value")
	colnames(blank)<-temp16
	
	rows<-c("Comparison")
	
	for(l in 1:(index-1)){rows[l+1]<-l}
	row.names(blank)<-rows
	
#	wilcoxcomps<-data.frame(blank)
	HTML.title("<bf>All pairwise comparisons - Asymptotic Mann-Whitney tests", HR=2, align="left")
	HTML(blank, classfirstline="second", align="left")
#	HTML(wilcoxcomps, classfirstline="second", align="left")
	HTML.title("<bf> ", HR=2, align="left")
	HTML.title("<bf>Why are there two different tests presented? The Behrens Fisher tests are the recommended approach, although we still need to independently verify these results.", HR=0, align="left")

} else if (statstest=="CompareToControl"&& leng >= 3)
{
	#Comparisons back to a control
	npmccontrol<-function(response,treatment, cont){
	dataset<-data.frame(var=response,class=treatment)
	summary(npmc(dataset, df=1, alpha=sig, control=contlevel), type="Steels")};
	steel <- npmccontrol(eval(parse(text = paste("statdata$", response))), eval(parse(text = paste("statdata$", treatment))), paste(contlevel))
	
	Treatment<-steel$'Data-structure'[2]
	Nobs<-steel$'Data-structure'[3]
	temp3<-cbind(Treatment, Nobs)
	header<-c(NA, " ")
	steeltable<-rbind(header, temp3)
	steeltable1<-data.frame(steeltable)
	temp4<-c("Treatment", "Nobs")
	colnames(steeltable1)<-temp4
	good2<-c("Treatment ID")
	for(i in 1:length(unique(eval(parse(text = paste("statdata$", treatment)))))){good2[i+1]<-i}
	row.names(steeltable1)<-good2

	#Steels Test results table
	HTMLbr()
	HTML.title("<bf>Steel's all comparisons back to one", HR=2, align="left")

	#Comparison<-steel$'Results of the multiple Steel-Test'[1]

	group2<-c("a")	
	for(i in 1:(length(unique(eval(parse(text = paste("statdata$", treatment)))))-1)){group2[i]<- paste(contlevel)}
	vs<-c("r")	
	for(i in 1:(length(unique(eval(parse(text = paste("statdata$", treatment)))))-1)){vs[i]<-"vs."}

	group1<-unique(eval(parse(text = paste("statdata$", treatment))))
	group1<-group1 [ group1!= contlevel]

	Pvalue<-format(round(steel$'Results of the multiple Steel-Test'[6],4), nsmall=4, scientific=FALSE)
	
	temp5<-cbind(group1,vs , group2, Pvalue)

	for (i in 1:(dim(temp5)[1])) 
	{
		if (temp5[i,4]<0.001) 
		{
			temp5[i,4]<-0.001
			temp5[i,4]<- paste("<",temp5[i,4])
		}
	}

	header<-c(NA, NA, NA,NA)
	temp7<-rbind(header, temp5)
	steeltable2<-data.frame(temp7)
	temp6<-c("Group", "vs.", "Group", "p-value")
	colnames(steeltable2)<-temp6
	good3<-c("Comparison")
	for(i in 2:length(unique(eval(parse(text = paste("statdata$", treatment)))))){good3[i]<-i-1}
	row.names(steeltable2)<-good3
	HTML(steeltable2, classfirstline="second", align="left")
}

#Text for conclusions, comments and references
textindex=1
HTMLbr()
HTML.title("<bf>Analysis conclusions", HR=2, align="left")

if (statstest=="MannWhitney" || (statstest=="AllComparisons" && leng == 2 ) || (statstest=="CompareToControl" && leng == 2 ) )
{
	dim <- length(unique(eval(parse(text = paste("statdata$", treatment)))))
	
	if (dim > 2)
	{
		if (temptab[1,3] <= (1-sig))
		{
			add<-paste(c("There is a statistically significant overall difference between the treatment groups at the "), sep="")
			add<-paste(add, 100*(1-sig), sep="")
			add<-paste(add, "% level of significance as the p-value is less than ", sep="")
			add<-paste(add, 1-sig, sep="")
			add<-paste( add, " (Kruskal-Wallis test)." , sep="")
			HTML.title(add, HR=0, align="left")
		} else if (temptab[1,3]>0.05)
		{
			add<-paste(c("The overall difference between the treatment groups is not statistically significant at the "), sep="")
			add<-paste(add, 100*(1-sig), sep="")
			add<-paste(add, "% level of significance as the p-value is greater than ", sep="")
			add<-paste(add, 1-sig,sep="")
			add<-paste( add, " (Kruskal-Wallis test)." , sep="")
			HTML.title(add, HR=0, align="left")
		}
		HTMLbr()
		HTML.title("<bf>Analysis description", HR=2, align="left")
		HTML.title("The overall treatment effect was assessed using the non-parametric Kruskal-Wallis test, see Kruskal and Wallis (1952, 1953).", HR=0, align="left")
		HTML.title("<bf> ", HR=2, align="left")
		HTML.title("Non-parametric tests should be used if the data is non-normally distributed, 
	the variability is different between treatment groups or the responses are not continuous and numerical. ", HR=0, align="left")
		HTMLbr()
		HTML.title("<bf>Statistical references", HR=2, align="left")
		HTML.title("Kruskal, WH and Wallis, WA (1952). Use of ranks in one criterion variance analysis. JASA, 47, 583-621.", HR=0, align="left")
		HTML.title("<bf> ", HR=2, align="left")
		HTML.title("Kruskal, WH and Wallis, WA (1953). Errata for Kruskal-Wallis (1952). JASA, 48, 907-911." , HR=0, align="left")
	}
	else if (dim == 2)
	{
		if (temptab[1,2] <= (1-sig))
		{
			add<-paste(c("The difference between the two treatment groups is statistically significant at the "), sep="")
			add<-paste(add, 100*(1-sig), sep="")
			add<-paste(add, "% level of significance as the p-value is less than ", sep="")
			add<-paste(add, 1-sig, sep="")
			#add<-paste(add, ", (p=", sep="")
			#add<-paste(add, temptab[1,2], sep="")
			add<-paste( add, " (Mann-Whitney test)." , sep="")
			HTML.title(add, HR=0, align="left")
		} else if (temptab[1,2]>0.05)
		{
			add<-paste(c("The difference between the two treatment groups is not statistically significant at the "), sep="")
			add<-paste(add, 100*(1-sig), sep="")
			add<-paste(add, "% level of significance as the p-value is greater than ", sep="")
			add<-paste(add, 1-sig, sep="")
			add<-paste( add, " (Mann-Whitney test)." , sep="")
			HTML.title(add, HR=0, align="left")
		}
		HTMLbr()	
		HTML.title("<bf>Analysis description", HR=2, align="left")

		HTML.title("The difference between the two treatments was assessed using the non-parametric Mann-Whitney test, see Wilcoxon (1945), Mann and Whitney (1947).", HR=0, align="left")
		HTML.title("<bf> ", HR=2, align="left")
		HTML.title("Non-parametric tests should be used if the data is non-normally distributed, 
	the variability is different between treatment groups or the responses are not continuous and numerical. ", HR=0, align="left")
		HTMLbr()
		HTML.title("<bf>Statistical references", HR=2, align="left")
		HTML.title("Mann, HB and Whitney, DR (1947). On a test of whether one of two random variables is stochastically larger than the other. Annals of Mathematical Statistics, 18, 50-60.", HR=0, align="left")
		HTML.title("<bf> ", HR=2, align="left")
		HTML.title("Wilcoxon, F (1945). Individual comparisons by ranking methods. Biometrics Bulletin, 1, 80-83.", HR=0, align="left")
	}
} else if (statstest=="AllComparisons" && leng >= 3)
{
	add<-c(" ")
	for(y in 2:int) 
	{
		if (allpairtab[y,4] <= (1-sig))
		{
			add<-paste(add, "The difference between groups ", sep="")
			add<-paste(add, allpairtab[y,1], sep="")
			add<-paste(add, " and ", sep="")
			add<-paste(add, allpairtab[y,3], sep="")
			add<-paste(add, " is statistically significant at the ", sep="")
			add<-paste(add, 100*(1-sig), sep="")
			add<-paste(add, "% level of significance as the p-value is less than ", sep="")
			add<-paste(add, 1-sig, sep="")
			add<-paste( add, " (Behrens Fisher test). " , sep="")
			textindex<-textindex+1
		}
	}

	textindex2<-1
	add2<-c(" ")
	for(q in 2:index)
	{
		if (blank2[q,1] <= (1-sig))
		{
			add2<-paste(add2, "The difference between groups ", sep="")
			add2<-paste(add2, blank[q,1], sep="")
			add2<-paste(add2, " and ", sep="")
			add2<-paste(add2, blank[q,3], sep="")
			add2<-paste(add2, " is statistically significant at the ", sep="")
			add2<-paste(add2, 100*(1-sig), sep="")
			add2<-paste(add2, "% level of significance as the p-value is less than ", sep="")
			add2<-paste(add2, 1-sig, sep="")
			add2<-paste(add2, " (Mann-Whitney test). " , sep="")
			textindex2<-textindex2+1
		}
	}
	if (textindex==1) 
	{
		HTML.title("<bf>None of the all pairwise Behrens Fisher tests were significant.", HR=0, align="left")
	} else if (textindex > 1) 
	{
		HTML.title(add, HR=0, align="left")
	}
	if (textindex2==1) 
	{
		HTML.title("<bf> ", HR=2, align="left")
		HTML.title("<bf>None of the all pairwise Mann-Whitney tests were significant.", HR=0, align="left")
	} else if (textindex2 > 1) 
	{
		HTML.title("<bf> ", HR=2, align="left")
		HTML.title(add2, HR=0, align="left")
	}
	HTMLbr()
	HTML.title("<bf>Analysis description", HR=2, align="left")
	HTML.title("All pairwise differences between the treatments were assessed using Behrens Fisher tests, see Munzel and Hothorn (2001) and Mann-Whitney tests, see Mann and Whitney (1947).", HR=0, align="left")
	HTML.title("<bf> ", HR=2, align="left")
	HTML.title("Non-parametric tests should be used if the data is non-normally distributed, 
the variability is different between treatment groups or the responses are not continuous and numerical. ", HR=0, align="left")
	HTMLbr()
	HTML.title("<bf>Statistical references", HR=2, align="left")
	HTML.title("Mann, HB and Whitney, DR (1947). On a test of whether one of two random variables is stochastically larger than the other. Annals of Mathematical Statistics, 18, 50-60.", HR=0, align="left")
	HTML.title("<bf> ", HR=2, align="left")
	HTML.title("Munzel, U and Hothorn, LA (2001). A unified approach to simultaneous rank test procedures in the unbalanced one-way layout. Biometrical Journal, 43(5) 553-569.", HR=0, align="left")
	HTML.title("<bf> ", HR=2, align="left")
	HTML.title("Wilcoxon, F (1945). Individual comparisons by ranking methods. Biometrics Bulletin, 1, 80-83.", HR=0, align="left")
} else if (statstest=="CompareToControl" && leng >= 3)
{
	add<-c(" ")
	for(y in 2:length(unique(eval(parse(text = paste("statdata$", treatment)))))) 
	{
		if (temp7[y,4] <= (1-sig))
		{
			add<-paste(add, "The difference between groups ", sep="")
			add<-paste(add, temp7[y,1], sep="")
			add<-paste(add, " and ", sep="")
			add<-paste(add, temp7[y,3], sep="")
			add<-paste(add, " is statistically significant at the ", sep="")
			add<-paste(add, 100*(1-sig), sep="")
			add<-paste(add, "% level of significance as the p-value is less than ", sep="")
			add<-paste(add, 1-sig, sep="")
			add<-paste(add, " (Steel's test). " , sep="")
			textindex<-textindex+1
		}
	}
	if (textindex==1) 
	{
		HTML.title("<bf>None of the Steel's all to one comparisons were significant.", HR=0, align="left")
	} else if (textindex > 1) 
	{
		HTML.title(add, HR=0, align="left")
	}
	HTMLbr()
	HTML.title("<bf>Analysis description", HR=2, align="left")

	HTML.title("The comparison of treatment groups back to a single control group was made using the non-parametric Steel's test, see Steel (1959).", HR=0, align="left")
	HTML.title("<bf> ", HR=2, align="left")
	HTML.title("Non-parametric tests should be used if the data is non-normally distributed, 
the variability is different between treatment groups or the responses are not continuous and numerical. ", HR=0, align="left")
	HTMLbr()	
	HTML.title("<bf>Statistical reference", HR=2, align="left")
	HTML.title("Steel, RGD (1959). A multiple comparison rank sum test: treatments versus control. Biometrics, 15, 560-572. ", HR=0, align="left")
}

HTMLbr()
HTML.title("<bf>R references", HR=2, align="left")

HTML.title("<bf> ", HR=2, align="left")
	HTML.title("<bf>   R Development Core Team (2008). R: A language and environment for
	statistical computing. R Foundation for Statistical Computing,
	Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.", HR=0, align="left")

HTML.title("<bf> ", HR=2, align="left")
	HTML.title("<bf>   Alan Genz, Frank Bretz, Torsten Hothorn with contributions by
	Tetsuhisa Miwa, Xuefei Mi, Friedrich Leisch and Fabian Scheipl
	(2008). mvtnorm: Multivariate Normal and t Distributions. R package
	version 0.9-0.", HR=0, align="left")

HTML.title("<bf> ", HR=2, align="left")
	HTML.title("<bf>     Joerg Helms and Ullrich Munzel (2008). NPMC: Nonparametric Multiple Comparisons. R package version
	1.0-7.
	", HR=0, align="left")

HTML.title("<bf> ", HR=2, align="left")
	HTML.title("<bf>    Lecoutre, Eric (2003). The R2HTML Package. R News, Vol 3. N. 3, Vienna, Austria.
	", HR=0, align="left")
